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Airfoil  Design  by  an  All-At-Once  Method  * * * § 

Ajit  Shenoy  t  Matthias  Heinkenschloss  *  Eugene  M.  Cliff 


Abstract 

The  all-at-once  approach  is  implemented  to  solve  an  optimum  airfoil  design  problem. 
The  airfoil  design  problem  is  formulated  as  a  constrained  optimization  problem  in  which  flow 
variables  and  design  variables  are  viewed  as  independent  and  the  coupling  steady  state  Eu¬ 
ler  equation  is  included  as  a  constraint,  along  with  geometry  and  other  constraints.  In  this 
formulation,  the  optimizer  computes  a  sequence  of  points  which  tend  toward  feasiblility  and 
optimality  at  the  same  time  (all-at-once).  This  decoupling  of  variables  typically  makes  the 
problem  less  nonlinear  and  can  lead  to  more  efficient  solutions.  In  this  paper  an  existing  op¬ 
timization  algorithm  is  combined  with  an  existing  flow  code.  The  problem  formulation,  its 
discretization,  and  the  underlying  solvers  are  described.  Implementation  issues  arc  presented 
and  numerical  results  are  given  which  indicate  that  the  cost  of  solving  the  design  problem  is 
approximately  six  times  the  cost  of  solving  a  single  analysis  problem. 

Key  words  Airfoil  design,  optimization,  computational  fluid  dynamics,  Euler  equa¬ 
tions,  nonlinear  programming,  optimal  design. 


1  Introduction 

Optimum  airfoil  design  is  an  active  area  of  research.  See,  e.g.,  the  recent  papers  [2],  [4],  [18], 
[19],  [20],  [22],  [23],  [26],  [27],  [33],  [34].  Abstractly,  the  optimum  airfoil  design  problem  can  be 
formulated  as  a  constrained  optimization  problem,  and  many  techniques  have  been  applied  to  its 
solution.  Most  of  the  recent  approaches,  in  fact  all  of  the  references  above,  combine  optimization 
and  optimal  control  techniques  with  computational  fluid  dynamics.  Several  issues  have  to  be  dealt 
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with  when  one  follows  this  path.  Among  these  are  the  formulation  of  the  optimization  problem, 
the  discretization  of  the  infinite  dimensional  airfoil  design  problem,  and  the  differentiation  of  ob¬ 
jective  function  and  constraints.  For  the  solution  of  the  airfoil  design  problem  all  issues  have  to  be 
dealt  with  simultaneously.  Great  care  must  be  taken  to  avoid,  or  at  least  to  control  inconsistencies 
and  to  develop  a  robust  and  efficient  solution  method.  While  these  techniques  have  already  been 
successfully  applied  to  airfoil  design  problems,  further  investigations  are  needed  to  improve  effi¬ 
ciency  and  robustness  of  the  solution  techniques  and  to  increase  the  set  of  (airfoil)  design  problems 
to  which  these  techniques  can  be  applied. 

In  this  paper  we  investigate  the  applicability  of  the  all-at-once  formulation  of  the  optimization 
problem  to  solve  an  airfoil  design  problem.  In  most  papers  on  airfoil  design,  the  flow  variables 
q  are  viewed  as  functions  of  the  design  parameters  w.  This  function  q(w)  is  implicitly  defined 
by  the  governing  equations  R(q1  w)  =  0,  in  our  case  the  steady  state  2-D  Euler  equations.  The 
optimization  formulation  describing  the  airfoil  design  problem  is  then  posed  in  the  design  variables 
w.  This  is  called  the  black-box  approach.  The  Euler  equations  are  not  visible  to  the  optimizer,  but 
hidden  by  eliminating  the  flow  variables,  i.e.,  by  expressing  the  flow  variables  q  as  functions  of 
the  designs  variables  w.  An  alternative  to  this  approach  is  the  all-at-once  formulation  in  which 
one  views  flow  variables  q  and  design  variables  w  as  independent  variables  in  the  optimization 
problem.  The  Euler  equations  coupling  these  two  are  included  into  the  optimization  formulation 
as  a  constraint  along  with  other  constraints  such  as  geometric  constraints,  drag  constraints,  etc.  The 
optimizer  is  now  responsible  for  computing  a  point  which  is  feasible  and  optimal  at  the  same  time, 
i.e.,  move  towards  feasibility  and  optimality  at  once,  rather  than  moving  along  the  manifold  of 
feasible  points  towards  optimality.  Comparisons  between  these  two  approaches  on  other  problems 
have  shown  that  the  all-at-once  approach  can  be  substantially  faster.  The  reason  is  that  viewing 
q  and  w  as  independent  variables,  allows  the  optimizer  to  violate  the  Euler  equations  during  the 
iterations.  These  are  only  required  to  be  satisfied  at  the  solution.  This  makes  the  optimization 
problem  less  nonlinear  and  often  results  in  fewer  iterations.  For  example,  Iollo  et  al.  [18]  report 
that  their  pseudo-time  stepping  implementation  of  the  all-at-once  approach  requires  only  three  to 
four  times  as  many  iterations  to  solve  the  design  problem  as  compared  to  the  effort  required  for 
the  solution  of  a  single  analysis  problem.  Our  results  indicate  a  factor  of  five  or  six.  However, 
our  optimization  approach  is  different  and  our  problem  includes  geometric  constraints.  It  is  also 
important  to  note  that  an  optimizer  implementing  the  all-at-once  approach  requires  roughly  the 
same  problem  information  as  an  optimizer  applied  to  the  black-box  approach,  except  that  the  all- 
at-once  approach  does  not  require  solutions  to  the  nonlinear  flow  equations.  We  give  a  more 
detailed  presentation  of  the  relations  in  the  next  section. 

Rather  than  formulating  the  airfoil  design  problem,  its  discretization,  and  a  solution  algorithm 
and  then  implement  all  components  from  scratch,  we  decided  to  build  upon  existing  codes.  In 
our  implementation  of  the  all-at-once  method  for  our  airfoil  design  problem  we  combine  the 
optimizer,  TRICE,  with  the  flow  code,  ErICA.  This  imposes  certain  limits  on  the  choice  of  problem 
formulation  and  discretization,  but  we  believe  this  to  be  a  realistic  approach.  As  we  have  indicated 
above,  various  issues  have  to  be  addressed  when  solving  the  airfoil  design  problem.  We  focus 
on  the  optimization  formulation.  Our  airfoil  parameterization  is  obtained  by  choosing  a  set  of 
basis  airfoils  and  computing  an  optimized  airfoil  as  a  linear  combination  of  those.  Moreover, 
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grid  generation  and  discretization  of  the  Euler  equation  was  done  to  limit  difficulties  arising  from 
nonsmoothness  and  inconsistencies.  Further  gains  in  efficiency  and  accuracy  can  be  achieved  by 
using  more  refined  discretization  techniques  and  improving  the  coupling  of  the  flow  solver  with 
the  optimizer.  This  was  beyond  the  scope  of  this  study  and  is  planned  for  future  investigations. 

This  paper  is  organized  as  follows:  In  Section  2  we  discuss  optimization  formulations  and  their 
relations.  This  section  also  provides  further  motivation  for  the  all-at-once  approach  and  reviews 
some  existing  optimization  approaches.  The  governing  equations  and  the  flow  code  ErICA  are 
discussed  in  Section  3.  The  design  problem  is  formulated  in  Section  4  and  Section  5  contains  a 
description  of  the  optimizer  TRICE.  Section  6  contains  some  implementation  issues  that  have  to 
be  resolved  when  combining  an  optimizer  with  a  flow  code  for  our  situation.  Section  7  presents 
some  numerical  results  and  contains  a  discussion  of  our  numerical  experiments  and  open  issues. 

2  Optimization  Problem 

There  are  several  ways  to  cast  the  design  problem  outlined  in  the  introduction  into  an  optimization 
problem.  Two  formulations  will  be  discussed  in  this  section.  The  main  purpose  of  this  section  is 
to  provide  a  background  for  the  discussion  of  our  approach  to  the  airfoil  design  problem  and  for 
a  comparison  with  other  approaches  in  the  literature.  In  this  section  we  proceed  as  follows:  First, 
we  present  the  two  formulations  and  their  relation  in  an  abstract  framework.  Then  we  discuss  the 
applicability  to  the  airfoil  design  problem. 

The  first  formulation  of  the  airfoil  design  problem  is  given  by 


min  J(q,w), 

(2.1) 

s.t.  R(q,  w)  =  0, 

(2.2) 

G(q ,  w)  <  0. 

(2.3) 

Here  q  represent  the  flow  variables  and  w  are  the  design  parameters.  The  constraint  function  R 
represents  the  Euler  equations.  The  inequality  constraints  (2.3)  represent  geometric  constraints, 
drag  constraints  and  the  like.  The  special  case  in  which  G  does  not  depend  on  the  flow  parameters 
q  deserves  attention.  In  this  section  we  assume  that  the  functions  J  :  IRnq  x  lRnw  — »  1R  ,  R  : 
IRnq  x  Nriv  — *  lRUq ,  and  G  :  JR"''  x  JRnw  — »  Kn,J  are  twice  continuously  differentiable  at  the  points 
under  consideration.  However,  we  note  that  for  the  formulation  and  execution  of  the  optimization 
algorithm  applied  to  our  airfoil  design  problem,  we  only  need  first  derivatives.  For  R,  G  we  denote 
differentiation  with  respect  to  a  variable  by  using  the  variable  as  a  subscript,  e.g.,  Rq(q,  w)  denotes 
the  partial  Jacobian  of  R  with  respect  to  q.  In  addition  to  the  differentiability  assumption,  we  make 
the  assumption  that  Rq(q,  w )  is  invertible  at  all  points  (q,  w)  under  consideration. 

Under  the  assumption  of  the  implicit  function  theorem  (see,  e.g.,  [24]),  the  constraint  (2.2) 
locally  defines  a  function  q  :  IRnw  — »  IRnq  as  the  solution  of 

R(q(w),w)  =  0.  (2.4) 

If  the  equation  (2.2)  has  a  unique  solution  q(w)  for  all  w  e  IRnw  under  consideration  (typically, 
(2.3)  represents  an  explicit  restriction  of  the  design  space  and  therefore  not  the  whole  IRnw  is 
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relevant),  then  we  can  eliminate  the  flow  variables  q  and  formulate  (2. 1)— (2.3)  in  the  following 
reduced  form: 


min  j(w)  =  J(q(w),w),  (2.5) 

s.t.  G(w)  =  G(q(w),w )  <  0.  (2.6) 

The  optimization  formulation  (2. 1)— (2.3)  corresponds  to  the  all-at-once  (AAO)  approach  [7], 
[11]  (also  called  the  simultaneous  analysis  and  design  (SAND)  approach  [3],  [26]).  The  optimiza¬ 
tion  formulation  (2.5),  (2.6)  corresponds  to  the  black-box  approach  [11]  (also  called  the  nested 
analysis  and  design  (NAND)  [3],  [26])  and  it  corresponds  to  the  multidiscipline  feasible  and  indi¬ 
vidual  discipline  feasible  approach  in  [7] . 1 

In  the  following,  we  present  optimality  conditions  for  (2. 1)— (2.3)  and  we  discuss  the  relation 
between  these  two  problems.  These  results  are  known  and  can  be  found  in  a  similar  form,  e.g.,  in 
[9],  [13].  Let 

L(q,  w,  A,  /j)  =  .J(q,  w)  +  A TR(q,  w)  +  RTG(q,  w)  (2.7) 

be  the  Lagrangian  corresponding  to  (2. 1)— (2.3).  If  a  constraint  qualification  is  met,  then  for  an 
optimal  point  (q,  w)  of  (2. 1)— (2.3)  there  exist  A,  //  such  that 

VqJ{q,w)  +  Rq{q,w)T\  +  Gq(q,w)T/J,  = 

VwJ(q,w)  +  Rw(q,w)T\  +  Gw(q,w)T/j  = 

R(q,w)  = 

G(q,w)  < 

E  > 

G(q,w)Tn  = 

If  G  does  not  depend  on  q,  then  the  first  equation  in  (2.8)  reduces  to 

V9  J(q,  w)  +  Rq(q,  w)T A  =  0. 

Equation  (2.9)  is  called  the  adjoint  equation  and,  if  G  does  not  depend  on  q,  defines  the  Lagrange 
multiplier  (or  the  co-state)  A.  With  A  given  by  (2.9),  the  term  VwJ(q1  w)  +  Rw(q1  w)TX  of  the 
second  equation  in  (2.8)  is  called  the  reduced  gradient. 

A  commonly  used  constraint  qualification  is  the  linear  independent  constraint  qualification 
(LICQ):  Let  G1  (q,  w)  denote  the  vector  of  functions  of  G(q1  w)  which  are  active  at  (q1  w).  Then 
LICQ  is  satisfied  if  the  gradient  of  the  component  functions  in  R(q.  w)  and  G'  (q.  w)  are  linearly 
independent.  If  G  does  not  depend  on  q  and  if  the  rows  of  G1  (w)  are  linear  independent,  which  is, 
e.g.,  the  case  if  G(w)  =  ±w,  then  our  assumption  that  Rq(q,  w)  is  invertible  implies  that  LICQ  is 
satisfied  [9],  [13]. 

The  second  order  necessary  [sufficient]  optimality  conditions  are  given  by  (2.8)  and 

(Sq)  H(q,w,  X,n)  (Sq^j  “  0  (2.10) 

\swJ  \swJ  [>j 

'Since  we  only  have  one  discipline  R(q,w)  =  0  there  is  no  distinction  between  the  multidiscipline  feasible  and 
individual  discipline  feasible  formulation  in  [7]. 


0, 

0, 

0, 

0, 

0, 

0. 


(2.8) 


(2.9) 
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for  all  sq,  sVJ  satisfying 


Rq(q,w)sq  +  Rw(q,w)sw  =  0,  (2.11) 

Gq(q,w)sq  +  Gw(q,w)sw  =  0.  (2.12) 

Unless  noted  otherwise,  G(q,  w )  would  usually  refer  to  the  set  of  active  constraints,  G1  (q,  w).  Here 
H(q,  w,  A,  //)  denotes  the  Hessian  of  the  Lagrangian 

H(q,  w,  A,  n)  =  V2{g,w)L(q,  w,  A,  /i). 

Points  satisfying  the  homogeneous  state  equation  (2. 11)  can  be  characterized  by  )  =  T(q,w)sw, 
where 

(2.13) 

With  this,  (2.10),  (2.11),  (2.12)  can  be  rewritten  as 

swT(q ,  w)TH(q ,  w,  A,  n)T(q,  w)sw  >  [>]  0  (2.14) 

for  all  sw  satisfying 

[~Gq(q,w)Rq(q,  w)~1Rw(q,  w)  +  Gw(q,w)]sw  =  0.  (2.15) 

The  matrix  T(q,  w)TH(q,  w,  A,  f^)T(q,  w)  is  called  the  reduced  Hessian.  The  term  on  the  left  hand 
side  of  (2.15)  can  either  be  computed  by  calculating  the  sensitivities  Rq(q1  w)~1Rw(q1  w)  or  using 
an  adjoint  approach.  If  we  define 

A  =  -Rq(q(w),  w)~TGq(q(w),  w)T ,  (2.16) 

then  the  left  hand  side  of  (2.15)  can  be  written  in  the  form  AT Rw(q(w)1  w)  +  GV](q(w)1w).  In 
particular  if  ng  is  smaller  than  nw  the  adjoint  equation  based  approach  seems  more  attractive  than 
the  sensitivity  equation  approach. 

It  is  known,  see,  e.g.,  [9],  [13],  that  derivatives  for  the  reduced  problem  (2.5),  (2.6)  are  related  to 
the  reduced  quantities  of  the  problem  (2. 1)— (2.3).  For  example,  the  gradient  V ./(  «.;)  of  the  reduced 
problem  is  equal  to  the  reduced  gradient  NwJ(q1  w)  +  Rw(q1  w)TX,  with  A  given  by  (2.9),  at  q  = 
q(w).  Moreover,  the  Hessian  H(w1  / 1 )  =  V2,L(wi  d)  °f  the  Lagrangian  L{w ,  /i)  =  J(w)+G(w)T /i 
of  the  reduced  problem  is  equal  to  the  reduced  Hessian  T ( q,  w)TH(q1  w1  A,  /i)T ( q1  w) .  Finally,  the 
Jacobian  Gw(w)  is  equal  to  the  matrix  on  the  left  hand  side  of  (2.15)  at  q  —  q(w). 

Comparisons  between  the  all-at-once  (SAND)  approach  and  the  black-box  (NAND,  discipline 
feasible)  approach  can  be  found  in,  e.g.,  [7],  [11],  [26].  The  all-at-once  approach  decouples 
state  and  design  variables.  An  optimizer  for  (2. 1)— (2.3)  can  use  this  decoupling  and  is  allowed 
to  violate  constraints  during  the  iteration.  This  can  result  in  substantial  gains  in  performance. 
Significant  reductions  in  solution  times  for  problems  related  to  the  one  considered  in  this  paper 
are  reported  in  [11],  [12].  However,  the  optimizer  must  achieve  feasibility  and  optimality  at  the 
same  time.  This  requires  carefully  designed  optimization  codes  to  maintain  robustness.  In  the 
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computational  experiments  in  [11]  the  (particular)  implementation  of  the  black-box  formulation 
always  performed  more  robustly  than  the  implementation  of  the  all-at-once  approach  for  a  one¬ 
dimensional  duct  design  problem.  Concerning  the  applicability  of  the  all-at-once  approach  and  the 
black-box  approach,  it  should  be  mentioned  that  most  of  the  quantities  needed  to  implement  the 
all-at-once  approach  also  have  to  be  provided  for  an  implementation  of  the  black-box  approach. 
This  is  indicated  by  the  relations  between  derivatives  for  the  reduced  problem  (2.5),  (2.6)  and  the 
reduced  quantities  of  the  problem  (2. 1)— (2.3)  summarized  above. 

In  the  context  of  airfoil  design  problems  both  formulations  (2. 1)— (2.3)  and  (2.5),  (2.6)  have 
been  used,  however,  currently  the  black-box  formulation  (2.5),  (2.6)  seems  to  be  dominant  [19], 
[20],  [23],  [27],  [33].  For  airfoil  design  problems  the  all-at-once  formulation  (2. 1)— (2.3)  is  con¬ 
sidered  in  [18],  [34].  In  both  cases  only  the  equality  constrained  problem  (2.1),  (2.2)  is  considered. 
The  optimization  methods  are  derived  from  the  optimality  system  for  (2.1),  (2.2).  In  [18]  the  opti¬ 
mality  system  is  solved  by  applying  a  few  pseudo-time  steps  to  the  Euler  equation  and  the  adjoint 
equation  to  improve  state  and  co-state  estimates  and  a  gradient-like  step  to  update  the  design  vari¬ 
ables.  The  optimization  methods  in  [34]  are  derived  from  the  application  of  Newton’s  method  to 
the  optimality  system;  these  are  particular  versions  of  sequential  programming  (SQP)  methods. 

We  use  the  all-at-once  formulation  (2. 1)— (2.3).  For  our  particular  design  problem  the  con¬ 
straints  G  are  simple  constraints  on  the  design  variables.  The  exact  problem  formulation  will  be 
introduced  in  the  subsequent  sections.  We  use  an  SQP  method  from  the  class  of  methods  described 
in  [9],  [15]  for  the  solution  of  the  all-at-once  formulation.  This  SQP  method  uses  an  interior  point 
strategy  to  handle  the  inequality  constraints  and  employs  a  trust-region  strategy  for  globalization 
of  convergence  and  to  enhance  robustness.  See  also  Section  5.  If  only  equality  constraints  are 
present,  then  the  Newton  based  methods  in  [34]  are  related  to  the  SQP  methods  in  [9],  [15].  Be¬ 
sides  the  capability  of  handling  inequalities  on  the  designs,  other  main  differences  are  that  the  SQP 
methods  in  [9],  [15]  use  a  trust-region  globalization  and,  in  addition  to  exact  second  derivatives, 
provide  quasi-Newton  approximations  to  the  full  and  reduced  Hessian  of  the  Lagrangian.  First 
and  second  order  convergence  results  are  proven  in  [9]  and  the  influence  of  inexact  derivatives  is 
analyzed  in  [14]. 

For  the  formulation  of  the  airfoil  design  problem  as  an  optimization  problem,  several  other 
issues  are  of  great  importance.  These  are  the  issues  of  discretization,  differentiability,  and  unique 
solvability  of  state  equations  and  linearized  state  equations.  We  give  a  more  detailed  description 
below.  For  general  airfoil  design  problems  comprehensive,  rigorous  treatments  of  these  issues 
are  still  missing.  In  the  case  of  a  one-dimensional  duct  design  problem,  which  is  related  to  the 
airfoil  design  problem,  such  a  comprehensive,  rigorous  treatment  can  be  found  in  [6].  It  is  shown 
in  [6]  that  an  understanding  of  these  issues  can  be  used  to  improve  robustness  and  efficiency  of 
the  optimization  code.  These  improvements  are  based  on  the  understanding  of  the  problem,  of  its 
discretization,  and  of  the  optimization  method.  They  are  achieved  with  very  little  programming 
effort  and  almost  no  additional  computing  effort  per  iteration. 

The  airfoil  design  problem  originally  is  an  infinite  dimensional  problem.  Therefore,  the  opti¬ 
mization  formulation  and  optimization  algorithm  have  to  be  combined  with  a  discretization  scheme. 
Various  approaches  are  possible.  Two  of  those  are  the  optimize-then-discretize  approach  in  which 
the  optimization  algorithm  is  formulated  in  the  infinite  dimensional  setting  and  then  discretization 
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are  applied  to  the  individual  steps,  and  the  discretize-then-optimize  approach  in  which  one  first 
discretizes  the  problem  and  then  applies  an  optimization  algorithm  to  the  discretized  problem.  The 
processes  of  discretization  and  optimization  are  usually  not  interchangeable  and  therefore  these 
two  approaches  are  different.  The  numerical  solution  of  an  infinite  dimensional  problem  requires 
a  careful  study  of  the  problem  at  hand.  Several  issues  have  to  be  kept  in  mind.  In  the  optimize- 
then-discretize  approach  the  derivatives  after  discretization  are  usually  not  the  derivatives  of  the 
discretized  functions.  Therefore  optimization  algorithms  have  to  cope  with  inexact  derivative  in¬ 
formation.  See  e.g.,  [5],  [33].  The  discretize-then-optimize  approach  often  neglects  the  fact  that 
the  infinite  dimensional  problem  structure  still  influences  the  finite  dimensional  problem.  If  this 
influence  is  not  incorporated  properly,  then  the  optimization  problem  typically  becomes  artificially 
ill-conditioned  and  one  observes  a  severe  degradation  in  performance  and  robustness  of  the  opti¬ 
mizer,  see  [6],  [28]. 

Derivatives  of  constraint  functions  and  solutions  q  to  the  state  equations  are  used  in  the  formu¬ 
lation  of  optimality  conditions  and  in  efficient  optimizers.  See,  for  example,  gradient  computations 
using  sensitivities  or  adjoint  equations.  For  problems  governed  by  the  Euler  equations,  differen¬ 
tiability  in  the  infinite  dimensional  context  is  problematic,  due  to  the  presence  of  shocks.  This 
might  be  different  for  the  discretized  Euler  equations.  If  smoothing  procedures  (e.g.,  introduction 
of  artificial  viscosity)  are  applied  in  discretization  schemes  for  the  Euler  equations,  the  resulting 
finite  dimensional  system  may  be  differentiable.  However,  since  the  discretization  schemes  used 
in  CFD  codes  are  very  complex,  ‘derivatives’  and  ‘adjoint  equations’  should  be  treated  with  care 
and  usually  must  be  understood  formally. 

It  is  also  important  to  keep  in  mind  that  the  formulations  (2. 1)— (2.3)  and  (2.5),  (2.6)  are  only 
equivalent  if  (2.2)  has  a  unique  solution  q(w )  for  all  w  e  IRnw  under  consideration.  If  R(q,  w)  =  0 
represents  the  (discretized)  Euler  equations,  this  assumption  seems  to  be  rather  strong  in  view  of 
the  non-uniqueness  result  presented  in  [21]  for  discretized  Euler  equations.  The  existence  and 
uniqueness  of  the  solution  q  of  R(q,  w)  =  0  for  given  w  is  often  also  related  to  the  existence  and 
uniqueness  of  the  solution  sq  of  the  linearized  state  equations  Rq(q ,  w)sq+Rw(q,  w)sw+R(q ,  w )  = 
0  for  given  ( q ,  w),sw. 

As  we  have  noted  before,  for  a  one-dimensional  duct  design  problem,  which  is  related  to  the 
airfoil  design  problem,  the  above  issues  have  been  rigorously  discussed  in  [6].  For  general  air¬ 
foil  design  problems  these  issues  are  subject  of  current  research.  In  our  approach  to  the  airfoil 
design  problem  we  parameterize  the  airfoil  using  linear  combinations  of  existing  airfoils.  This 
can  be  viewed  as  a  reduced  basis  approach  (see  e.g.,  [32,  Sec. 7. 3])  leading  to  a  low  dimensional 
(nw  =  4)  design  space.  Our  grid  generation  scheme  leads  to  grids  which  depend  smoothly  on  the 
design  parameters.  Our  application  programs  are  based  on  the  package  ErICA  for  the  simulation  of 
flows  over  airfoils  governed  by  the  Euler  equations.  Among  the  discretization  schemes  available 
in  that  package,  we  use  the  schemes  with  better  smoothness  properties.  We  use  the  discretize- 
then-optimize  approach.  Since  we  have  a  low  dimensional  design  space  and  a  rather  simple  grid 
generation  scheme,  we  believe  this  is  sensible.  However,  given  our  experiences  in  [6],  we  believe 
this  approach  has  to  be  rethought  if  more  complex  discretization  schemes  are  used.  More  details 
on  the  discretization  schemes  are  provided  in  Sections  3  and  6. 
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3  Analysis  Problem,  Discretization,  and  Flow  Code 

In  this  section  we  discuss  the  analysis  problem  underlying  our  design  problem  and  its  discretiza¬ 
tion.  The  analysis  problem  is  the  flow  q  around  the  airfoil  governed  by  the  steady  state  Euler 
equations  for  a  perfect  gas.  We  also  outline  the  flow  code  ErICA  used  for  the  solution  of  the  analy¬ 
sis  problem.  Although  our  optimization  formulation  is  based  on  the  all-at-once  approach  and  our 
optimizer  never  needs  to  solve  the  Euler  flow  equations,  we  will  extract  several  subtasks  from  the 
flow  code  ErICA.  The  presentation  of  the  ErICA  code  will  help  to  describe  these  tasks. 

The  unsteady  Euler  equations  for  a  perfect  gas,  written  in  integral  conservation  law  form  is 
given  by 

f  QdS+  f  F-hds  =  0  (3.1) 

at  J n  Jsn 


where,  in  Cartesian  coordinates, 


F  —  T j  -Y  Qk 


pu 2  +  p 
puv 
(. ph0)u 


pv2  +  p 

(i pho)v 


with  velocity  components  u,  v,  density  p,  total  energy  per  unit  mass  eQ  —  e  +  (u2  +  v 2) /2,  with  e 
being  the  internal  energy  per  unit  mass  and  pressure  p,  which  for  a  perfect  gas  may  be  expressed 
by  the  relation,  p  =  ( 7  —  1  )pe.  The  total  enthalpy  per  unit  mass  is  given  by  ha  =  a2/ (7  —  1)  + 
u2 / 2  +  v2 / ‘2  =  ec  +p/p,  where  a  =  \Jjp]~p  is  the  sonic  velocity.  Here,  Q  represents  the  conserved 
variables  with  q  =  [p  u  v  p]T  denoting  the  primitive  variables,  and,  T  and  Q  represent  the  inviscid 
fluxes.  The  problem  domain  is  denoted  by  El  and  SQ  represents  the  boundary  of  the  domain.  For  a 
detailed  treatment  of  the  Euler  equations  refer  [10]. 


3.1  Discretization  of  the  Euler  Equations 

For  a  given  airfoil  configuration,  the  shape  of  which  is  represented  as  a  function  of  the  design 
variables  w,  the  analysis  problem  corresponds  to  the  solution  of  the  Euler  equations  of  flow.  The 
flow  is  simulated  numerically  using  the  solver  ErICA  (EuleR  Inviscid  Code  for  Aerodynamics) 
which  was  developed  by  Narducci  [25]. 

Computational  simulations  were  performed  on  a  C-type  grid,  which  is  wrapped  around  the 
airfoil.  We  only  sketch  the  grid  generation  to  fix  some  notation.  For  the  omitted  details  we  refer  to 
[29,  Sec.  4.2.2].  The  grid  is  generated  algebraically,  by  the  following  procedure:  We  first  distribute 
points  on  bottom  boundary,  corresponding  to  the  airfoil  surface  and  the  trailing  edge  wake,  and 
on  the  top  boundary  of  the  computational  grid,  corresponding  to  the  far-field  boundary.  Once  the 
boundaries  are  defined,  we  connect  corresponding  pairs  of  points  on  the  top  and  bottom  boundaries 
using  straight  lines.  Grid  cells  and  nodes  are  numbered  by  (j,  k),  where  j  refers  to  the  horizontal 
position  and  k  refers  to  the  vertical  position  in  the  grid.  Indices  with  k  =  1  refer  to  nodes  or  cells 
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on  or  at  the  airfoil,  respectively.  A  typical  201  x  53  grid  is  shown  in  Figures  3.1  and  3.2  with  121 
points  on  the  airfoil  surface.  In  practical  CFD  codes  more  sophisticated  grid  generation  schemes 
are  used.  Eventually,  such  grid  generations  have  to  be  incorporated.  However,  in  a  first  attempt  to 
apply  the  all-at-once  methodology  to  airfoil  design,  we  preferred  this  simple  grid  because  of  the 
relative  ease  of  generating  the  grid  and  because  of  its  guaranteed  smooth  dependence  on  the  design 
parameters  (airfoil). 


Figure  3.1:  The  201  x  53  grid. 


Figure  3.2:  Close-up  View  of  the  201  x  53  grid. 


Given  the  grid,  the  ErICA  code  [25]  is  used  for  the  solution  of  the  steady  state  Euler  equations. 
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In  our  version  of  ErICA  a  finite  volume  discretization  using  an  upwind  scheme  with  Van  Leer  Flux 
Vector  Splitting  is  applied  to  compute  the  residuals.  A  MUSCL  (Monotone  Upstream-centered 
Scheme  for  Conservation  Laws)  differencing  approach  is  used  to  interpolate  the  values  of  the  state 
variables  q  from  the  cell  centers  to  the  cell  faces.  In  order  to  fully  capture  the  shock,  we  use 
third  order  interpolation  of  the  fluxes.  We  use  Van  Albada’s  limiter  to  suppress  oscillations  in 
the  flow  solution.  While  other  discretization  schemes  are  available  in  ErICA,  we  selected  these 
because  of  their  better  smoothness  properties.  See  the  discussion  at  the  end  of  Section  2.  A 
pseudo-time  marching  scheme  is  used  to  compute  solution  to  the  steady  state  Euler  equations. 
We  use  the  boundary  conditions  in  [16,  Ch.  19]  on  the  airfoil  surface  and  the  far- fie  Id  boundary 
conditions  proposed  in  [30].  We  describe  the  main  features  of  the  flow  solver  which  are  needed  in 
the  subsequent  discussion.  For  more  details  we  refer  to  [29]. 

As  mentioned  above,  we  use  a  cell-centered,  finite  volume  formulation  to  rewrite  the  governing 
equations  (3.1).  For  the  (j,  /c)th  grid-cell  the  (semi-discretized)  residual  in  terms  of  the  primitive 
variables  is  given  by 

do 

+  Rjk(q,  w)  =  0,  (3.2) 

where  S]k  is  the  area  of  the  ( j ,  k) th  subdomain,  M  =  AQ  is  the  Jacobian  of  the  mapping  between 

the  conserved  and  the  primitive  variables,  and  Rjk  =  Y.Sides(F  ■  n)As  is  the  residual,  where  F  is 
the  inviscid  flux  and  As  is  the  length  of  the  side.  Summation  is  done  over  all  sides  of  cell  (j,  k ). 
Requiring  (3.2)  for  all  cells  yields 


do 

SM-A  +  R{q,w)  =  0.  (3.3) 

The  residual  is  computed  as 

”0  =  [C  ■  *>H-1/2+ [(i? '  fl>H+1/2+ [T  ■  ftH*-1/2+IV  ■  flHt+1/2  ■  <3-4> 

where  [F-n]j± i/2  correspond  to  the  inviscid  flux,  F ■ h ,  across  the  vertical  cell  faces,  and  [F-h]k±i/2 
corresponds  to  the  flux  across  the  horizontal  cell  faces,  respectively.  For  a  given  cell  face,  we  have, 


pU 

puU  +  hxp 

pU  V  +  flyP 

(. pho)U 


where  U (=  hxu  +  hyv)  is  the  velocity  normal  to  the  cell  face,  and  hx  and  hy  are  the  Cartesian 
components  of  the  normal  to  the  cell  face.  As  noted  above,  the  residual  is  computed  using  an 
upwind  scheme,  with  Van  Leer  Flux  Vector  Splitting,  with  third  order  interpolation  via  MUSCL 
differencing  to  interpolate  the  values  of  the  state  variables,  q,  from  the  cell  centers  to  the  cell  faces, 
and  with  Van  Albada’s  limiter  to  suppress  oscillatory  behavior  in  the  flow  solution. 

We  use  the  surface  boundary  conditions  described  in  [16,  Ch.  19]  and  [8]  (flow  tangency  stip¬ 
ulation  and  requirements  that  the  normal  momentum  be  zero  and  entropy  be  conserved,  curvature 
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corrections  described  in  [8]  are  neglected),  along  with  the  far-held  boundary  conditions  proposed 
in  [30].  These  boundary  conditions  involve  the  creation  of  ghost  cells. 

The  steady  state  solution  corresponding  to  the  semi-discrete  equations  (3.3)  is  computed  using 
a  pseudo-time  marching  scheme  applied  to  (3.3).  We  consider  the  implicit  scheme 


A  q  =  -Rn, 


(3.5) 


where  A q  =  qn+ 1  —  qn;  qn  =  q(nAt').  Here  ^  denotes  an  approximation  of  the  Jacobian  F  R. 
See  below.  The  scheme  (3.5)  may  be  regarded  as  a  simplified  implicit  Euler  scheme  since, 


SM 


A  q 
A t 


-K 


>n+l 


-  R" 


-Rn 


A  q. 


The  equation  (3.5)  is  ‘solved’  by  applying  one  step  on  an  Alternating  Direction  Implicit  (ADI) 
scheme.  The  resulting  time-marching  scheme  is  not  accurate  in  time,  but  this  is  not  required  since 
we  are  only  interested  in  steady  state  solutions. 

For  a  moment,  suppose  that  {§^R)n  =  {§^R)n  and  that  A  =  Fp.  Then  the  equation  of  (3.5) 
corresponding  to  cell  (j.  k )  is  given  by 

{^M + (M  ■  ">H— 1/2 + M  •  ")H+V2 

+  [(i  ■  + 1*-4  ■  ftH*+1/2) }  =  -R(&-  (3-6> 

See  (3.2),  (3.4).  Instead  of  using  A  =  A.F,  we  make  two  simplifications  to  derive  A.  These  in¬ 
crease  the  efficiency  with  which  one  step  of  the  pseudo-time  marching  scheme  can  be  performed. 
The  first  simplification  is  as  follows.  In  the  residual  computation,  flux  terms  like  [(F  ■  n)  As]J+1/2 
are  calculated  using  Van  Leer  Flux  Vector  Splitting  and  MUSCL  differencing  with  cubic  interpo¬ 
lation  of  the  values  of  the  state  variables  q  from  the  cell  centers  to  the  cell  faces.  ErICA  analyt¬ 
ically  computes  the  Jacobians  of  the  flux  terms  obtained  using  linear  instead  of  cubic  interpola¬ 
tion.  The  second  simplification  in  computing  A  is  made  by  partly  suppressing  the  influence  of  the 
ghost-cells.  Emphasizing  the  influence  of  the  boundary  conditions,  the  residual  can  be  written  as 
R(q1  w)  =  R(q1  qg(q),w),  where  qg  are  the  values  of  the  flow  variables  on  the  ghost  cells.  The 
boundary  conditions  are  used  to  express  these  as  functions  of  the  flow  variables  q  in  the  interior 
and  of  w:  qg  =  qg(q,  w).  Thus,  the  derivative  of  the  residual  is  of  the  form 

dR(q,  w)  =  dR(q,qg,w )  dR(q,qg,w )  dqg(q,w ) 

dq  dq  dqg  dq 

In  ErICA  the  approximation 

dR(q,w )  ^  dR(q,qg,w ) 
dq  dq 


(3.8) 
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is  used  to  obtain  A. 

Let  A  fa  Ap  denote  the  approximate  flux-Jacobians  derived  using  the  two  simplifications  out¬ 
lined  above.  If  we  set  A  =  (jA.  then  A  is  a  block  pentadiagonal  matrix.  We  split  A  =  Ak  +  Aj, 
where  Aj  corresponds  to  the  terms  [(A  •  n)Asl  +  [(A  •  ft)  As]  in  (3.6)  and  Ak  corre- 

sponds  to  the  terms  [(A  •  n)As]  +  [(A  •  n)Asj  in  (3.6).  The  subscript  j  in  Aj  indi¬ 
cates  that  the  matrix  includes  information  of  A  along  constant  j -lines  (vertical  grid-lines).  Simi¬ 
larly,  Ak  includes  information  of  A  along  constant  /c-lines  (horizontal  grid-lines).  We  also  define 
T  =  A-tSM.  Now,  (3.5)  can  be  written  as 

[T  +  Ank+  A]]  A q  =  -Rn.  (3.9) 

We  do  not  solve  (3.9),  but  approximately  factor 

T  +  Ank+A^[T  +  Ak]n  T-1  [T  +  Aj]n 

and  solve 

-I  r  -I 

T  +  Ak\  T-^T  +  Aj]  A q  =  -Rn  (3.10) 

The  step  A q  is  computed  by  solving 

[Tr+-4‘i"f„9l/2  =  (3.11) 

[T  +  Aj]  A  q  =  TAq1/2. 

The  new  flow  iterate  is 

qn+l  =  qn  +  Aq.  (3.12) 

The  two  simplifications  in  the  flux-Jacobians  leading  to  A  ^  guarantee  that  (after  sym¬ 
metric  permutation)  the  matrices  on  the  left  hand  sides  of  (3.11)  are  block  tridiagonal.  Thus,  each 
subproblem  in  (3.11)  requires  a  block  tridiagonal  matrix  inversion,  which  involves  a  block  LU 
factorization  and  a  block  matrix  solve;  the  latter  consists  of  forward  and  backward  substitutions. 
For  more  details  about  the  scheme,  we  refer  to  Hirsch  [17,  16]  or  Shenoy  [29].  An  outline  of  the 
ErICA  algorithm  for  the  solution  of  the  governing  Euler  flow  equations,  R(q,  w)  =  0  for  given  w, 
is  given  in  Algorithm  3.1. 

3.2  Airfoil  Shape  Parameterization 

We  use  an  airfoil  shape  parameterization,  formulated  in  [31].  See  also  [32,  Sec.  10].  The  airfoil 
geometry  is  represented  as  the  weighted  combination  of  six  shape  functions 

6 

y{x/c)  =  J2wiyi(x/c)-  (3.13) 

i— 1 

Four  of  the  shape  functions  are  pre-existing  airfoils,  namely,  NACA  2412,  NACA  64i-412,  NACA 
652-415  and  NACA  642-A215.  The  shape  functions  y\-y\  may  be  referenced  from  [1],  The  two 
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Algorithm  3.1  (ErICA) 

1  Given  Wi. 

1.1  Generate  C-grid,  including  ghost  cells. 

1.2  Compute  direction  cosines  and  lengths  for  each  cell  face  and  the  areas  of  each  cell. 

2  Given  qn.  Compute  residual: 

2. 1  Impose  Boundary  Conditions. 

2.2  Compute  Rjk  =  £sides(.F  '  fi)As. 

2.3  Compute  ||i?||. 

2.4  If  \\R\\  <  tol,  then  output  the  result  and  stop;  otherwise  goto  3. 

3  Euler  Implicit  Time  Integration. 

3.1  k  =  constant  lines: 

Solve  [T  +  Aky  Aq1/2  =  -Rn. 

3.2  j  =  constant  lines: 

Solve  | T  +  AjJ  A q  =  TAq1/2. 

3.2  Update  State:  qn+1  =  qn  +  A qn.  Set  n  =  n  +  1  and  goto  2. 


additional  shape  functions  2/5,  y6  are  used  to  impose  certain  geometric  closure  conditions  at  the 
trailing  edge  of  the  airfoil.  These  shapes  are  given  by 


3/5 


y& 


x/c,  on  upper  surface, 

0,  on  lower  surface, 

0,  on  upper  surface, 
— x/c,  on  lower  surface. 


Since  these  functions  are  used  to  close  the  airfoil  at  the  trailing  edge,  the  weights  and  w6  are 
fixed  in  terms  of  w\-wt±.  We  require  that 


2/us(l)  0; 


which  yield  the  following  relations 

®5  =  -blus(l)^1+l/2us  (1)^2  +  1/3US  (1)^3 +  t/4us(l)^4], 

=  [yhs  (l)wi  +  y2ls  (l)ro2  +  Vsu  (1)1373  +  2/4,.  (1)1374], 

where  subscripts  us  and  Is  refer  to  the  upper  and  lower  surfaces,  respectively.  An  efficient  approach 
to  implement  the  above  is  to  close  each  airfoil  j/i— 2/4  individually  to  obtain  yi-y4,  and  use  these  as 


14 


A.  R.  SHENOY,  M.  HEINKENSCHLOSS,  AND  E.  M.  CLIFF 


our  design  bases.  We  have 

4 

y(x/c)  =  Yj^iViix/c). 

i— 1 


(3.14) 


4  The  Design  Problem 

Given  {w*},  the  flow  q  around  the  airfoil  is  governed  by  the  Euler  equations  for  a  perfect  gas.  The 
design  problem  is  formulated  as  follows: 


max  C'l(q.  vj)  (4.1) 

such  that 

R(q,m)  =  0,  (4.2) 

CD(q,w)  <  CDm^,  (4.3) 

•S'min  <  S(w)  <  Smax,  (4.4) 

e{xu)  > 

^min-  (4.5) 


where  w  =  {zui},  and  q  =  [p  u  v  p]T  denote  the  primitive  variables  of  flow,  with  the  usual 
notation.  Equation  (4.2)  refers  to  the  discretized  steady  state  Euler  equations  of  flow.  The  drag, 
Cd,  in  this  case,  is  the  wave  drag.  The  lower  limit  on  the  area,  .S',  is  imposed  so  that  the  airfoil 
does  not  become  too  thin,  a  requirement  for  structural  integrity.  The  upper  limit  is  imposed  to  avoid 
thick,  unrealistic  airfoils.  Equation  (4.5)  represents  a  bound  on  the  trailing  edge  angle  imposed  to 
avoid  situations  in  which  the  upper  surface  can  go  below  the  lower  surface.  See  below.  The  free- 
stream  conditions  are  based  on  Mach  number  M  =  0.75  flow  at  angle  of  attack  a  =  0. 

The  state  equation  (4.2)  was  discussed  in  Section  3.  We  briefly  describe  the  computation  of  the 
aerodynamic  forces  used  to  compute  Cl(<7,  w)  and  Cn(q.  w ),  and  the  trailing  edge  condition  (4.5). 
These  are  fairly  standard,  but  are  included  for  completeness. 

The  aerodynamic  forces  are  computed  by  numerically  integrating  the  pressure  over  the  surface 
of  the  airfoil.  The  normalized  forces  normal  and  tangential  to  the  airfoil  chord  line  are  respectively 
given  by, 


2  ^ 

Cn  ~~  —  T7  ^  Pjk{Xj+l,l  ~ 

Poo  j=jo 

2  ^ 

CT  =  77  53  Pjk(yj+l,l  -  Vj, l), 

Poo  Voo  j  jo 


where  (j,  k ),  k  =  1,  j  =  j0, . . . ,  jj ,  denote  the  grid  points  on  the  airfoil  surface.  We  compute  the 
lift  and  drag  forces  respectively  as, 


Cl  =  Cn  cos  a  —  CT  sin  a, 
Cd  =  Cn  sin  a  +  Ct  cos  a. 


(4.6) 
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The  dependence  of  the  lift  and  drag  coefficients  on  the  state  q  and  the  design  w  can  be  determined 
from  (4.6). 

The  area  of  the  airfoil  (for  unit  chordlength)  is  given  by 

S  =  [  Vnsdx  -  [  yisdx. 

Jo  Jo 

Here  yus,  y\s  represent  the  upper  and  the  lower  surface,  respectively.  Representing  the  airfoil  in 
terms  of  the  basic  (closed)  airfoils  (3.14),  we  have 


S  = 


'uj-iS-i. 


where  S)  correspond  to  the  areas  of  the  individual  airfoils.  The  areas  of  the  individual  airfoils  can 
be  computed  at  the  beginning  of  the  design  cycle.  For  given  w  the  area  is  then  simply  computed 
as  the  weighted  sum  of  the  areas  of  the  given  (closed)  airfoils. 

Our  parametric  representation  of  the  airfoil  (3.14)  allows  for  situations  where  the  upper  surface 
can  go  below  the  lower  surface  of  the  airfoil.  Such  situations  were  actually  encountered  in  our 
preliminary  attempts  at  optimization.  Hence,  we  need  to  impose  an  additional  constraint  to  prevent 
such  physically  incompatible  configurations  to  arise.  This  is  done  by  constraining  the  trailing  edge 
angle  of  the  airfoil.  The  trailing  edge  angle  is  given  by  tan_1(y1's(l))  —  tan_1(y(IS(l)),  which  is 
approximated  by 

dxE  =  2/(s(l)  -  f/usl1)- 

Using  the  (approximate)  trailing  edge  angles  5^Ei  of  the  four  basic  airfoils,  this  can  be  written  as 


4 

<5t E  =  ^  TXJjdTEi- 
i— 1 


(4.7) 


It  was  found  sufficient  to  impose  the  requirement  (4.5)  to  ensure  that  the  upper  surface  does  not 
go  below  the  lower  surface.  Note  that  while  the  above  approximation  of  the  trailing  edge  angle  is 
fairly  crude,  it  does  yield  a  constraint  that  is  easy  to  compute  and  achieves  the  desired  effect. 


5  Optimization  Algorithm 

The  optimization  algorithm  used  for  our  computation  is  a  version  of  the  trust-region  interior-point 
SQP  methods  called  TRICE  developed  in  [9],  [14],  [15]  for  solving 


min  J(q,w), 

(5.1) 

s.t.  R(q,  w)  =  0, 

(5.2) 

WW  <  w  <  wmax. 

(5.3) 
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Clearly,  (5.1)— (5.3)  is  a  particular  case  of  (2. 1)— (2.3).  In  this  section  we  give  a  brief  description  of 
the  algorithm.  We  leave  out  many  technical  details  and  focus  on  how  the  algorithm  interfaces  with 
the  flow  solver.  For  more  details  on  the  algorithm  and  its  convergence  we  refer  to  the  papers  [9], 
[14].  An  important  aspect  of  the  TRICE  implementation  [15]  is  that  application  specific  subtasks 
in  the  optimization  are  separated  from  the  optimizer  TRICE  and  can  be  provided  by  the  user.  In 
our  case  this  allows  us  to  provide  approximate  solutions  to  linearized  state  equations  and  adjoint 
equations  computed  by  a  modification  of  the  flow  code  ErICA.  See  also  Section  6.  Since  (5.3) 
corresponds  to  (2.3)  with  a  G  independent  of  q,  the  Lagrange  multiplier  A  is  determined  by  (2.9). 
We  use  the  equation  (2.9)  to  define  A  =  A (q,w). 

If  we  define  a  diagonal  scaling  matrix  D(q,  w)  €  HE1"'  %n',!  with  diagonal  elements 


,\  (ww  ~w)f  if  (T(q,w)TVJ(q,w)).<  0, 
[D(q,w)  ..  =  <  i  ) 

{w-wmin)f  if  (T(q,  w)TVJ(q,  «>)).>  0, 

then  the  first  order  optimality  conditions  (2.8)  can  be  equivalently  written  as 


R(q,w)  =  0, 

D(q,  w)2T(q,  w)TVJ(q,  w)  =  0, 


(5.4) 


(5.5) 


and  wmin  <w<  wmax. 

The  algorithms  in  [9],  [14],  [15]  generate  a  sequence  of  iterates  (qklwk),  where  wk  is  strictly 
feasible  with  respect  to  the  bounds,  i.e.,  <  wk  <  wmax  (hence  the  term  interior-point 

method).  The  algorithms  can  be  motivated  by  applying  Newton’s  method  to  the  system  of  nonlin¬ 
ear  equations  (5.5)  where  the  w  component  is  kept  strictly  feasible  with  respect  to  the  bounds,  i.e., 
wmin  <  w  <  wmax.  The  step  s  =  (sq,  sw )  is  decomposed  in  to  a  quasi-normal  step  sn  and  a  tan¬ 
gential  step  st .  The  role  of  the  quasi-normal  step  sn  is  to  move  towards  feasibility.  It  is  of  the  form 
sn  =  (s“,  0).  The  ^-component  sq  is  related  to  the  Newton  step  applied  to  solve  R{q ,  wk)  =  0, 
for  given  wk.  The  role  of  the  tangential  step  is  to  move  towards  optimality.  It  is  of  the  form 
st  =  T(qk,  wk)sw  =  ( ~Rq(qk ,  wk)~lRw{(ik,  wk)sw,  sw),  where  T(qk,  wk)  is  the  representation  of 
the  null-space  of  the  linearized  state  equation  defined  in  (2.13).  The  u—component  sw  of  st  is 
related  to  a  quasi-Newton  step  for  the  reduced  problem  (2.5),  (2.6). 

The  matrix  D(q,w)  is  in  general  not  differentiable,  but  this  nondifferentiability  is  benign  and 
does  not  interfere  with  the  fast  convergence  of  Newton’s  method.  A  linearization  of  (5.5)  around 
qk,wk  gives 


( Rq)k&q  T  (-Riu) kSw 


~Rk, 

—  DkTkVJk. 


(5.6) 

(5.7) 


Here  we  have  used  the  subscript  k  to  denote  evaluation  of  functions  at  qk,wk.  In  (5.7),  0  de¬ 
notes  the  nw  x  nq  matrix  with  zero  entries,  £(q,  w.  A)  =  J(q,  w)  +  A Tc(q,  w),  Vfqw^(q,  w,  A)  = 
[J(q,  w)  +  A Tc(q,  u>)],  and  E(q,  w)  is  the  ]Rn^xn^  diagonal  matrix 

(E(q,w])..  =  \  {T(q,w)TNJ(q,w)). 
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replacing  the  in  general  not  existing  term  [d^u^  D2(q1  w)]T(q1  w)T'VJ(q1  w). 

Since  the  solution  of  the  linearized  state  equation  (5.6)  can  be  written  as 

S  —  S  T  Tk  j 

where  sn  =  (—(Rq^Rk,  0)T  and  Tk  is  given  by  (2.13). 

By  using  (5.8)  we  can  rewrite  the  linear  system  (5.6)-(5.7)  as 

s  —  s  T  Tksw, 

{pkTk  S7xxlkTkDk  +  Ekj  Dk  1sw  =  —DkTk  (N2q  w^£ksn  +  V  J/c), 


(5.9) 

(5.10) 


The  Newton-like  step  now  is  the  solution  of  (5.9),  (5.10)  with  Dk  replaced  by  Dk,  where  Dk  is 
defined  by  (5.4)  with  T£VJk  replaced  by  Tj [Vj^  +  VJ*].  This  change  of  the  diagonal 

scaling  matrix  is  based  on  the  form  of  the  right  hand  side  of  (5.10). 

One  can  that  if  ( qk,wk )  is  close  to  a  nondegenerate  minimizer  (q*,w*)  which  satisfies  the 
second  order  sufficient  optimality  conditions,  the  matrix  on  the  left  hand  side  of  (5.10)  is  posi¬ 
tive  definite.  Therefore,  (5.10)  can  also  be  interpreted  as  the  optimality  condition  of  a  quadratic 
program  in  sw.  To  globalize  the  convergence  and  to  enhance  robustness  of  the  algorithm,  a  trust- 
region  globalization  is  added.  Let  Ak  be  the  trust  radius  at  iteration  k.  The  ^-component  of  sn  is 
computed  by  approximately  solving 


minimize  4||(JR(?)/c(sn)(?  +  i4||2 
subject  to  ||(sn)q||  <  A*. 


(5.11) 


Given  sn,  the  step  in  w  is  computed  by  approximately  solving 

minimize  (t£  (. Hksnk  +  V  Jk))T  sw  +  (Tj HkTk  +  EkD AT2) 

subjectto  ||T)fc1s«;||  <  5k. 


(5.12) 


Of  course,  we  also  have  to  require  that  the  new  iterate  is  in  the  interior  of  the  box  constraints. 
To  ensure  that  wk  +  sw  is  strictly  feasible  with  respect  to  the  box  constraints  we  choose  ok  G 
[<j,  1),  <7  G  (0, 1),  and  compute  sw  with  ak(wm in  —  wk)  <  sw  <  ak(wm in  —  wk).  The  quadratic 
minimization  problems  (5.11)  and  (5.12)  only  needs  to  be  solved  approximately.  For  example,  an 
approximate  solution  of  (5.1 1)  is  given  by 


(S’1)?  =  ~t[(Rq)k 


(5.13) 


where  t  =  1  if  ||[(JR<?)/c]“1i4||  <  Ak  and  t  =  A/c/||[(i?g)fe]“17?fe||  <  Ak  otherwise. 

An  approximate  solution  of  (5.12)  can  be  computed  using  a  modified  conjugate  gradient  method. 
If  the  reduced  Hessian  T[ HkTk  is  approximated  by  a  quasi-Newton  update,  then  the  solution  of 
(5.12)  is  relatively  inexpensive.  The  cost  of  computing  the  ‘reduced’  gradient  T^{Hks\  +  V-4) 
dominates  the  cost  of  solving  (5.12). 

The  main  steps  of  the  trust-region  interior-point  SQP  scheme  are  outlined  in  algorithm  5.1. 
More  general  versions  of  these  algorithms  are  possible.  See  [9],  [14],  [15]. 
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Algorithm  5.1  (TRICE) 

1  Given  H0  and  A0. 

2  For  k  =  0, 1,  2, ...  do 

2.1  Compute  si  by  approximately  solving  (5.11). 

2.2  Compute  Tk ( Hksl  +  VJ*). 

2.3  Computes^  with(jA;(tymin  — <  sw  <  crk(wmaK—wk)  by  approximately  solving 
(5.12). 

2.4  Compute  s  =  sn  +  s*  =  sn  +  Tksw. 

2.5  Compute  A  by  solving  (2.9)  with  q  =  qk  +  sq,  w  =  wk  +  sw. 

2.6  Update  the  trust  region  radius  Ak  and  decide  if  qk  +  sq,  wk  +  sw  can  be  accepted 
as  the  new  iterate. 

2.7  If  s  is  rejected  set  qk+1  =  qk,  wk+i  =  wk,  and  Xk+i  =  Xk. 

Otherwise,  s  is  accepted;  set  qk+ 1  =  qk  +  sq,  wk+ 1  =  wk  +  sw  and  A/j+i  =  A. 

2.8  If  exact  second  order  information  is  not  used,  update  the  (reduced)  Hessian  ap¬ 
proximation. 


We  briefly  sketch  the  information  that  the  SQP  algorithm  5.1  requires  from  the  application 
programs.  A  detailed  presentation  is  given  in  [15].  If  (5.13)  is  used,  then  step  2.1  requires  the 
solution  of  a  linearized  state  equation.  The  computation  in  step  2.2  involves  the  solution  of  an 
adjoint  equation,  see  the  definition  (2.13)  of  T(q1  w).  If  a  quasi-Newton  approximation  is  used 
to  replace  the  the  reduced  Hessian  T^HkTk,  then  step  2.3  can  be  implemented  very  efficiently 
using  a  modified  conjugate  gradient  method.  The  application  of  Tk  in  step  2.4  requires  the  solu¬ 
tion  of  another  linearized  state  equation.  See  (2.13).  Computations  involving  the  solution  of  one 
linearized  state  equation  and  one  adjoint  equation  may  be  needed  in  step  2.8.  This  depends  on  the 
update  used.  This  and  other  algorithms  are  implemented  in  [15].  Global  and  local  convergence 
results  are  proven  in  [9].  The  influence  of  inexact  derivatives  is  analyzed  in  [14].  The  latter  aspect 
is  important  in  our  application  since  we  use  a  pseudo-time  marching  (iterative)  scheme  to  com¬ 
pute  approximate  solutions  to  linearized  state  equations  and  to  the  adjoint  equations.  Moreover, 
additional  approximations,  outlined  in  Section  3,  are  applied  in  these  computations  as  well. 

A  final  remark  on  the  handling  of  inequality  constraints  is  in  order.  In  our  case  the  design 
space  is  small  nw  —  4  and  other  approaches  such  as  projection  methods  or  active  set  methods  can 
likely  be  used  with  similar  performance  to  handle  the  inequality  constraints  (5.3).  However,  wing 
designs  in  industrial  settings  may  involve  up  to  500  design  parameters  [4].  In  this  case  an  interior 
point  approach  promises  to  be  superior. 
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6  Numerical  Implementation 

6.1  Handling  the  Inequality  Constraints  and  Reformulation  of  the  Opti¬ 
mization  Problem 


The  current  version  of  TRICE  only  solves  problems  of  the  form  (5. 1)— (5.3).  Hence  we  need  to 
recast  the  design  problem  (4.1)-(4.4)  into  the  form  (5. 1)— (5.3).  This  is  done  by  handling  (4.3)  as  a 
“soft”  constraint  using  a  penalty  term  and  by  transforming  the  design  parameters. 

Instead  of  including  the  drag  constraint  (4.3)  we  add  a  penalty  term  P(g{CD/CDmax  —  1))  to 
the  objective,  where 


P(z) 


0  z  <  0, 
z2  z  >  0. 


(6.1) 


Here  g  is  a  (scalar)  penalty  constant  which  can  be  used  to  increase  emphasis  on  the  drag  violation. 
The  addition  of  P(z)  to  the  objective  function  J  has  the  desired  effect  of  penalizing  the  objective 
when  the  drag  constraint  is  violated.  Note  that  this  is  a  “soft”  constraint,  in  the  sense  that  the 
optimizer  will  allow  the  drag  constraint  (4.3)  to  be  violated  as  long  as  the  penalty  term  added  is 
not  too  large.  This  can  be  addressed  to  some  extent  by  controlling  the  penalty  constant  g. 

The  remaining  constraints  can  be  addressed  using  a  mapping  between  the  control  variable,  w, 
and  the  design  weights,  w.  Rather  than  use  the  design  weights  as  our  control  variables,  we  use  the 
area  of  the  airfoil  and  its  trailing  edge  angle  as  control  variables,  as  shown  below.  This  enables  us 
to  address  the  issue  of  ensuring  that  the  lower  bound  on  the  area  and  the  trailing  edge  angle  remain 
strictly  enforced,  by  making  use  of  the  fact  that  we  can  place  bounds  on  the  control  variables.  We 
use,  as  our  control  variables, 


w  =  W  { 


2  •  (wi  —  VJ3) 
2  •  ("072  —  W4) 

2  •  S 


(6.2) 


0.5  •  <5t e 


where 


S  = 


Q  .  ’ 
^mtn 


5tk  — 


I JTE 
7  5 
Omin 


and  the  scalar  factor  W  has  been  added  in  order  to  be  able  to  experiment  with  the  scaling.  Note  that 
this  results  in  a  control  space  that  is  simply  a  result  of  a  linear  combination  of  the  design  weights, 
as  should  be  evident  from  equations  (4.7)  and  (4.7).  Now,  enforcing  bounds 


2W  <  w3  <  2W  • 

in 


strictly  enforces  the  bounds  on  the  area  of  the  airfoil  (4.4).  Similarly,  the  bound, 


0.5 W  <  W4 


imposes  the  constraint  (4.5).  The  mapping  (6.2)  can  be  used  to  translate  between  w  and  zu.  Note 
that  the  above  mapping  was  chosen  so  as  to  yield  a  low  condition  number  for  the  transformation 
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matrix  providing  the  mapping  between  the  weights  and  the  controls.  This  was  done  with  the 
philosophy  that  a  change  in  the  controls  should  produce  roughly  the  same  amount  of  change  in  the 
weights.  This  choice  did  yield  improved  performance  in  the  optimization  algorithm. 

The  original  design  problem  (4.1)-(4.4)  is  now  recast  as 

min9^  J(q,w)  =  ~CL(q,w)  +  P(g(CD(q,w)/CDmax  -  1)) 
s.t.  R(q,  w)  =  0, 

^min  A  W  E:  ®maxj 

which  is  of  the  form  (5.1)-(5.3). 

6.2  Solution  of  the  Linearized  State  and  Adjoint  Equations 

To  solve  the  design  problem  using  Algorithm  5.1  described  above,  we  need  to  be  able  to  do  the 
following: 

•  Provide  an  update  sq  in  the  state  variable,  given  an  update  in  the  control  variable  sw  by 
solving  the  linearized  state  constraint. 

•  Solve  the  adjoint  equation  at  a  given  point. 

These  tasks  can  be  performed  using  the  same  problem-solving  structure  applied  in  the  flow  solver 
ErICA.  These  tasks  can  be  easily  extracted  from  the  flow  code  and  only  relatively  few  changes  are 
needed. 

The  modification  of  the  ErICA  code  to  compute  approximate  solutions  to  the  linearized  state 
equation  is  shown  in  Algorithm  6.1.  The  scheme  outlined  in  Algorithm  6.1  solves  an  approxima¬ 
tion 

(a?  +Ak+  Ag)  Sq  +  ^(g,  w)sw  +  R(q ,  w )  =  R(sqi  sw,  q,w)  =  0  (6.3) 

of  the  linearized  Euler  equations  using  a  pseudo-time  marching  scheme  analogous  to  the  one  ap¬ 
plied  in  ErICA.  Here  q,  w  and  sw  are  given  and  an  approximate  solution  sq  has  to  be  computed. 
Note  that  §^R(q:  w)  is  replaced  by  Aj  +  Ak+Ag.  While  in  the  residual  computation,  flux  terms  are 
calculated  using  Van  Leer  Flux  Vector  Splitting  and  MUSCL  differencing  with  cubic  interpolation 
of  the  values  of  the  state  variables  q  from  the  cell  centers  to  the  cell  faces,  only  linear  interpola¬ 
tion  is  used  to  calculate  approximate  Jacobians,  see  Section  3.  This  leads  to  the  matrices  A].  Ak. 
However,  boundary  conditions  are  included  in  the  residual  computations,  cf.  (3.7).  This  is  reflected 
above  by  the  matrix  Ag.  The  equation  (6.3)  is  solved  by  driving  an  unsteady  form  of  the  linearized 
Euler  equations 

Os  ~  ~  ~  OR  — 

SM~0f  =  —(Aj  +  Ak  +  Ag)sq  -  —  {q,w)sw  -  R(q,w )  =  -R{sq,  sw,  q,w)  =  0  (6.4) 

towards  steady  state.  The  factor  SM  is  added  to  the  transient  term  in  order  to  make  the  above 
equation  consistent  with  the  discretized  Euler  equations  ( cf.  (3.3)).  The  pseudo-time  marching 
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scheme  used  is  identical  to  the  one  used  in  marching  the  nonlinear  Euler  equations,  described  above 
in  equations  (3. 10)— (3. 12),  with  the  nonlinear  residual,  R  replaced  by  the  linearized  residual,  R. 
We  can  view  this  algorithm  as  simply  an  iterative  method  for  solving  the  linearized  state  equation. 
Note  that  a  relaxation  factor  [5  is  used  to  update  the  solution  in  the  iterative  process.  Our  numerical 
experiments  showed  that  using  (5  =  1.25  yielded  improved  convergence  rates.  Also,  note  that 
there  is  an  external  loop  in  the  iterative  process  monitoring  the  residual.  This  is  in  order  to  ensure 
that  the  residual  does  not  diverge.  If  the  norm  of  the  residual  is  greater  than  some  predetermined 
value,  jRmaX’  then  the  iterative  process  is  restarted  with  a  reduced  time  step  At.  This  is  necessitated 
by  the  fact  that  the  Jacobians  can  be  ill-conditioned  if  the  solution  is  far  from  feasible,  resulting 
in  a  divergent  iteration.  Reducing  the  time  step  had  the  effect  of  alleviating  the  ill-conditioning. 
Using  -Rmax  =  1 2 1 1  JR0 1 1  seemed  adequate  for  our  purposes.  A  factor  of  12  is  used  because  in  some 
instances,  the  iterative  process  initially  increased  the  residual  but  managed  to  recover.  Clearly, 
these  rules  are  somewhat  ad-hoc  and  more  sophisticated  techniques  could  have  been  applied  to 
increase  efficiency.  Since  we  are  concerned  with  more  fundamental  issues  arising  in  the  all-at-once 
approach,  optimizing  performance  is  beyond  the  scope  of  this  study.  Note  that  Algorithm  6.1  only 
involves  one  Jacobian  evaluation  and  one  nonlinear  residual  evaluation.  The  Jacobian  undergoes 
one  block  LU  factorization  and  the  iterative  loop  only  involves  block  matrix  solves,  and  evaluation 
of  the  linearized  residual,  which  simply  requires  relatively  cheap  block  matrix  multiplications  and 
additions. 

Similarly,  for  the  solution  of  the  approximate  adjoint  equation 

(Aj  +  Ak  +  Ag)  A  +  ENq  J(q,  w )  =  A (q,  w)  =  0  (6.5) 

consider  a  “pseudo”  time  dependent  adjoint  equation, 

SMT  —  =  —  +  Ak  +  Agj  A  +  V qJ(zk)^  =  —A, 

where  J  is  the  objective.  As  in  the  solution  of  the  linearized  state  equation,  we  replace  -j^R  by 

Aj  +  Ak  +  Ag.  We  use  the  approximate  factorization  algorithm  used  in  ErICA  to  iterate  this 
equation  in  time,  until  the  residual  of  the  adjoint  equation  is  sufficiently  small,  ideally  A  =  0.  The 
procedure  is  as  follows.  We  have, 


| T  +  Aj  +  Ak J  AA  —  — An 

where  Aj  and  Ak  are  Jacobian  terms.  The  matrix  of  the  left  is  factored  approximately  according  to 
spatial  directions 

[T  +  Aj] T  T~t  [T  +  Ak] T  AA  =  -  An 
This  system  is  solved  using  the  sequence 

[7^  fV  =  “A"  (6.6) 

[T  +  Ak]  AA  =  TT  AA1/2 
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Algorithm  6.1  (Linearized  State  Equation  Solver) 

1  Given  q,  w,  sw,  tol. 

1.1  Generate  grid. 

1.2  Compute:  Aj(q,  w ),  Ak(q,  w),  Ag(q ,  w),  Rw(q ,  w). 

1.3  Compute:  R°  =  R(q1  w )  +  Rw(q,  w)sw. 

1.4  Set:  n  =  0,  sqn  =  0. 

2  LU  Decomposition. 

2.1  Compute:  LkUk  =  [T  +  Afej . 

2.2  Compute:  Lj  U3  =  |^T  +  Ajj . 

3  Euler  Implicit  Time  Integration. 

3.1  k  =  constant  lines: 

Solve  LkUkAsqi/2  =  -Rn. 

3.2  j  =  constant  lines: 

Solve  LjUjAsq  =  T  ■  A sqi/2. 

3.2  Update  sq:  s”+1  =  s”  +  A  sq.  Set  n  —  n+1  and  goto  4. 

4  Compute  Linearized  Residual 

4.1  Compute:  Rn  =  \A3  +  Afe  +  Ag j  sgn  +  R° 

4.2  Compute:  ||^n||. 

If  ||  Rn  1 1  <  tol,  set  sq  =  s^.  Return. 

Else,  if  ||^n||  <  Rmax,  goto  3. 

Else,  restart.  Set:  At  =  0.5Af,  n  =  0,  sqn  =  0,  and  goto  2. 


An+1  =  An  +  (3A\ 


The  above  iteration  is  performed  until  the  residual  of  the  adjoint  equation,  A,  is  reduced  to  zero. 
The  adjoint  computation  is  outlined  in  Algorithm  6.2.  Note,  that  while  computing  the  residual  of 
the  adjoint  equation,  and  the  linearized  state  equation,  we  include  the  terms  Ag  corresponding  to 
the  boundary  conditions. 
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Since  we  are  interested  in  the  solution  of  the  steady  state  adjoint  equation  (6.5),  we  could  have 
just  as  well  reversed  the  sequence  above,  i.e.,  solve 

[T  +  ^fAAv2  =  -A”  (6?) 

[t  +  Aj]  aa  =  TT  aa1/2 

and  set  An+1  =  An  +  [3  AX.  This  would  give  us  an  algorithm  which  is  exactly  the  same  as  that 
used  in  the  linearized  state  algorithm,  except  the  matrices  would  be  transposed.  However  even 
though  the  pseudo-time  marching  is  just  an  iterative  scheme  for  solving  (6.5),  we  preferred  to  use 
the  transpose  of  the  pseudo-time  process  (6.4)  to  calculate  the  adjoints.  Numerical  experiments 
showed  that  (6.6)  had  a  slightly  superior  convergence  behavior  than  (6.7). 

Once  again  an  outer  loop  monitors  divergence  of  the  residual.  We  use  Amax  =  12||A°||.  Nu¬ 
merical  experiments  showed  that  the  computation  of  the  adjoint  was  more  susceptible  to  producing 
divergent  results,  and  hence  care  has  to  be  taken  in  choosing  the  value  of  the  relaxation  factor  (3. 
We  choose  /3  =  min(1.25, 1  —  0.12  log(10 1| A® ||)),  which  has  the  desired  effect  of  underrelaxing 
when  the  solution  is  crude,  in  order  to  reduce  the  possibility  of  divergence,  and  overrelaxing  when 
the  solution  is  refined  in  order  to  increase  speed  of  convergence.  Also,  note  that  we  do  not  start  with 
A  =  0.  Rather,  we  start  from  the  previously  computed  estimate  of  the  adjoint  variable.  Our  exper¬ 
iments  showed  that  this  yielded  significant  savings  in  terms  of  the  number  of  iterations.  However, 
if  the  iteration  proves  to  be  divergent,  then  we  reset  A  to  zero.  As  we  have  noted  already  for  the 
linearized  state  solver,  these  rules  are  somewhat  ad-hoc  and  more  sophisticated  techniques  could 
have  been  applied  to  increase  efficiency.  This  will  be  done  in  future  studies.  As  with  the  procedure 
for  the  linearized  state  equation,  this  iteration  only  involves  a  single  Jacobian  evaluation. 


7  Numerical  Results  and  Discussion 

This  section  reports  on  some  of  numerical  experiments  conducted  using  the  TRICE  interior-point 
trust-region  SQP  optimization  algorithm  (see  Section  5)  coupled  with  the  modification  of  the  Er- 
ICA  flow  code  (see  Sections  3  and  5)  to  solve  the  airfoil  design  problem  stated  in  Sections  4  and  6. 
The  presentation  of  results  is  followed  by  a  discussion  of  observed  difficulties,  possible  remedies, 
and  further  research  issues. 

All  reported  computations  were  performed  on  C-type  grids  with  the  following  dimensions: 

1.  51  x  14  grid  with  31  points  on  the  airfoil  surface, 

2.  101  x  27  grid  with  61  points  on  the  airfoil  surface, 

3.  151  x  40  grid  with  91  points  on  the  airfoil  surface, 

4.  201  x  53  grid  with  121  points  on  the  airfoil  surface, 

5.  301  x  79  grid  with  181  points  on  the  airfoil  surface. 
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Algorithm  6.2  (Adjoint  Equation  Solver) 

1  Given  q,  w,  tol. 

1.1  Generate  grid. 

1.2  Compute:  Aj(q,  w ),  Ak(q,  w),  Ag(q ,  w),  Rw(q ,  w). 

1.3  Compute:  A0  =  \/qJ(q,w). 

1.4  Set:  n  =  0,  Xn  =  Aprev. 

(Aprev  is  the  Lagrange  multiplier  estimate  at  the  previous  iteration) 

2  LU  Decomposition. 

2.1  Compute:  LkUk  =  [T  +  A*]. 

2.2  Compute:  Lj U3  =  [T  +  Ay] . 

3  Compute  Adjoint  Residual 

3.1  Compute:  An  =  \A3  +  Ak  +  Ag j  An  +  A0 

3.2  Compute:  ||An||. 

If  ||  An  1 1  <  tol,  set  A  =  An.  Return. 

Else,  if  ||An||  >  Amax,  restart.  Set:  At  =  0.5Af,  n  =  0,  Xn  =  0,  goto  2. 

4  Euler  Implicit  Time  Integration. 

4.1  j  =  constant  lines: 

Solve  LjUjAX1/2  =  — An. 

4.2  k  =  constant  lines: 

Solve  LkUkA\  =  T  ■  AAi/2- 

4.3  Update  Adjoint:  An+1  =  Xn  +  AA.  Set  n  =  n  +  1  and  goto  3. 


Before  running  the  optimization,  we  used  ErICA  for  simulation  and  performed  a  grid  conver¬ 
gence  study.  This  was  done  for  the  (closed)  NACA  2412  airfoil  and  the  grids  specified  above.  The 
results  are  shown  in  Figure  7.1.  They  show  that  the  analysis  code  ErICA  overestimates  the  drag 
Cd  and  underestimates  the  lift  Cl  on  the  coarser  grids. 

The  grid  convergence  study  motivated  the  following  procedure  for  solving  the  design  problem. 
We  solve  the  design  problem  on  a  sequence  of  grids,  using  the  optimal  solution  (design  parameters 
w  and  interpolation  of  corresponding  flows  q)  of  the  coarse  grid  as  the  initial  value  on  the  next 
finer  grid.  On  the  coarse  grid  we  compute  starting  values  as  follows:  The  initial  design  parameters 
w  =  (1,  0,  0,  0)  correspond  to  the  (closed)  NACA  2412  airfoil  and  the  initial  flow  q  was  the  corre¬ 
sponding  flow  computed  by  ErICA  on  the  coarse  grid,  i. e. ,  we  start  with  a  point  satisfying  the  state 
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equation  R(q,  w )  =  0.  We  also  start  with  a  loose  bound  CDmax  on  the  drag,  which  is  tightened 
gradually  as  we  refine  the  grid.  The  values  are  reported  in  Table  7.1.  We  use  g  =  10  in  the  penalty 
term  for  the  drag  and  Smin  =  0.075  and  ,5max  =  0.15. 

We  stop  the  optimization  on  the  current  grid  if  the  norm  of  the  reduced  gradient  T(q,  w)  rVJ(q,  w ) 
is  reduced  to  a  tolerance  of  10-2  and  if  \\  R(q-  ) |  <  10-5.  It  was  observed  that  the  optimizer 
struggled  to  diminish  the  reduced  gradient  much  below  this.  For  the  51  x  14  grid,  a  scale  factor 
of  W  =  1  was  used,  which  produced  results  that  had  reduced  the  norm  of  the  reduced  gradient  to 
about  1.5  x  10-2,  at  which  point  the  iteration  was  terminated  due  to  lack  of  progress.  The  algo¬ 
rithm  was  implemented  so  that  if  the  optimization  stalled,  i.  e. ,  the  algorithm  terminated  because  it 
failed  to  make  further  progress,  we  restart  the  optimization  with  a  feasible  solution  for  the  given 
configuration.  Several  restarts  were  required  for  the  51x14  grid.  In  particular,  the  optimizer  strug¬ 
gled  if  the  constraint  residual  ||i?(<7,  w)  ||  became  too  large,  i.e.,  if  the  iterates  move  too  far  away 
from  feasibility.  The  observed  difficulties  in  the  optimization  are  likely  due  to  the  inaccuracies  of 
residual  Jacobians  used  in  the  linearized  state  equations  and  the  adjoint  equations.  We  will  discuss 
this  in  more  detail  below.  For  the  other  grids  we  used  a  scale  factor  of  W  =  10  (cf.  (6.2)).  The 
optimization  converged  with  a  single  restart  in  these  cases.  The  results  are  shown  in  Table  7.1. 
Note  that  we  did  not  run  the  optimizer  for  the  final  grid,  as  we  had  a  nearly  converged  solution. 
The  area  of  the  optimized  airfoil  is  at  the  lower  bound,  as  one  might  expect. 


Figure  7.1:  Grid  Convergence  Study  for  NACA  2412 

Figure  7.2  shows  the  final  results  obtained  for  the  301  x  79  grid.  Note  that  the  shock  has 
moved  downstream.  The  optimization  improves  the  performance  of  the  airfoil.  Compared  to  the 
baseline  NACA  2412  airfoil,  which  had  a  lift  coefficient,  Cl  =  0.4577  with  Cd  =  0.01241,  the 
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Table  7.1:  Numerical  Results  for  Airfoil  Design  Problem 


Grid  Size 

CDmax 

Data 

W1 

TU2 

W3 

TU4 

Cl 

Cd 

S 

51  x  14 

0.04 

Initial 

1.0000 

0.0000 

0.0000 

0.0000 

0.4132 

0.0469 

0.0823 

Final 

0.2538 

0.1793 

-0.1369 

0.5793 

0.5263 

0.0403 

0.0770 

101  x  27 

0.02 

Initial 

0.2538 

0.1793 

-0.1369 

0.5793 

0.5722 

0.0251 

0.0770 

Final 

0.2812 

0.1049 

-0.0760 

0.5314 

0.5227 

0.0201 

0.0750 

151  x  40 

0.014 

Initial 

0.2812 

0.1049 

-0.0760 

0.5314 

0.5266 

0.0159 

0.0750 

Final 

0.3059 

0.0593 

-0.0313 

0.5006 

0.5037 

0.0142 

0.0750 

201  x  53 

0.012 

Initial 

0.3059 

0.0593 

-0.0313 

0.5006 

0.5044 

0.0126 

0.0750 

Final 

0.3153 

0.0429 

-0.0154 

0.4893 

0.4955 

0.0120 

0.0750 

301  x  79 

0.010 

0.3153 

0.0429 

-0.0154 

0.4893 

0.4959 

0.0103 

0.0750 

final  design  has  Cl  =  0.4959  and  Cn  =  0.01028,  which  means  the  lift  coefficient  has  increased 
by  approximately  8.5%,  while  the  drag  coefficient  has  been  reduced  by  17%. 

Table  7.2  gives  an  account  of  the  computational  effort  required  at  each  step  to  produce  a  “con¬ 
verged”  solution.  However,  one  should  keep  in  mind  that  we  did  not  optimize  the  implementation 
for  efficiency.  Improvements  in  performance  can  be  achieved  and  we  outline  a  few  possible  en¬ 
hancements  below.  Here  the  number  of  “successful”  iterations  is  the  number  of  iterations  in  which 
the  Algorithm  5.1  accepts  the  step  s  and  uses  qk+1  =  qk  +  sq,  wk+\  =  wk  +  sw  as  the  new  iterate 
(see  step  2.7),  whereas  the  total  iterations  also  includes  the  unsuccessful  iterations,  i.e.,  those  in 
which  the  step  is  rejected  and  the  next  iterate  is  set  to  be  qk+ 1  =  qk,  wk+ 1  =  wk.  It  should  be 
noted  that  most  of  the  computational  effort  (in  terms  of  number  of  iterations)  occurs  at  the  coarsest 
level,  where  the  computations  are  fairly  cheap.  Though  there  is  room  for  improvement  the  TRICE 
algorithm  is  relatively  efficient  in  finding  the  given  solutions.  Consider,  for  example,  the  compu¬ 
tational  effort  required  for  the  51  x  14  grid.  We  require  3786  residual  evaluations,  5004  Jacobian 
evaluations,  10008  block  LU  factorizations  and  454948  block  matrix  solves.  Compare  this  to  the 
effort  required  to  obtain  a  single  analysis  solution:  We  require  approximately  1000  pseudo-time 
integration  steps  to  produce  a  converged  solution  which  requires  1000  residual  evaluations,  1000 
Jacobian  evaluations,  2000  block  LU  factorizations  and  2000  block  matrix  solves.  Discounting 
the  discrepancy  in  the  number  of  solves,  the  computational  effort  required  by  TRICE  is  roughly 
equal  to  the  effort  required  to  perform  5-6  flow  analyses,  which  is  very  cheap.  It  should  be  noted 
that  though  we  require  a  large  number  of  block  matrix  solves,  this  just  involves  forward  and  back¬ 
ward  substitutions  which  is  fairly  cheap.  As  we  have  indicated  earlier,  this  paper  concerned  with 
the  feasibility  of  the  all-at-once  approach  for  airfoil  design.  Computational  efficiency  was  not  a 
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a.  Surface  Pressure  Distribution 


b.  Initial  airfoil  (dashed)  and  optimal  airfoil  (solid) 


Figure  7.2:  Results  obtained  using  TRICE  for  Airfoil  Design  Problem 


prime  concern  for  the  interface  of  the  ErICA  flow  subroutines  with  the  TRICE  optimizer.  Sev¬ 
eral  improvements  can  be  made.  For  example,  we  recompute  the  Jacobian  information  every  time 
that  we  require  to  solve  the  linearized  state  equation  or  the  adjoint  equation.  Since  the  Jacobian 
will  only  change  if  the  iterate  (q,w)  changes,  a  more  efficient  implementation  would  require  only 
one  Jacobian  evaluation  for  each  iteration  of  Algorithm  5.1.  Moreover,  our  pseudo-time- stepping 
scheme  for  solving  the  linearized  state  and  the  adjoint  equation  can  be  improved  which  would  lead 
to  fewer  pseudo-time-steps  and  fewer  LU  solves.  Similar  instances  of  implementing  the  code  in  a 
more  efficient  manner  should  yield  significant  savings  from  the  elimination  of  redundant  compu¬ 
tations. 

Earlier,  we  have  pointed  out  that  the  optimizer  struggled  to  achieve  the  desired  tolerances.  We 
now  discuss  possible  reasons  for  this  and  also  possible  remedies.  We  distinguish  among  three 
groups.  The  first  group  is  related  to  the  information  that  is  provided  to  the  optimizer,  the  second 
group  is  related  to  the  optimization  formulation,  and  the  third  group  is  related  to  the  airfoil  design 
problem  and  its  discretization. 

The  TRICE  optimizer  requires  from  the  user  the  solution  of  linearized  state  equation  and  ad- 
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Table  7.2:  Computational  History  for  Airfoil  Design  Problem  using  TRICE 


Grid 

Total 

Successful 

Number 

#  Residual 

#  Jacobian 

#  LU 

#  Solves 

Size 

Iterations 

Iterations 

of  Restarts 

Evaluations 

Evaluations 

Factorizations 

51  x  14 

388* 

211 

12 

3786 

5004 

10008 

454948 

101  x  27 

44 

26 

1 

313 

437 

874 

73588 

151  x  40 

27 

11 

1 

506 

518 

1036 

81953 

201  x  53 

19 

5 

1 

386 

425 

850 

137768 

Iteration  was  stopped  with  a  norm  of  the  reduced  gradient  of  1.5  •  10  2  due  to  lack  of  progress. 


joint  equation.  We  extract  this  information  from  the  ErICA  code.  However,  we  use  simplifications 
in  the  computation  of  the  constrain  residuals.  While  the  residual  R(q,  w )  is  evaluated  using  Van 
Leer  Flux  Vector  Splitting  and  MUSCL  differencing  with  cubic  interpolation  of  the  values  of  the 
state  variables  q  from  the  cell  centers  to  the  cell  faces,  the  ‘Jacobians’  of  the  flux  terms  are  ob¬ 
tained  using  linear  instead  of  cubic  interpolation.  Thus  we  only  use  approximate  Jacobians.  These 
approximations  become  better  as  the  grid  is  refined,  but  on  a  given  grid  a  certain  error  level  cannot 
be  removed.  This  explains  our  earlier  observation  that  the  optimizer  had  more  problems  finding 
a  solution  relative  to  the  given  tolerances  on  the  coarse  grid  than  on  the  fine  grid.  The  Jacobians 
used  are  only  asymptotically  correct  and  the  discrepancies  between  true  Jacobians  and  Jacobians 
used  become  smaller  as  the  mesh  is  refined.  On  coarse  grids,  we  try  to  oversolve  the  problem. 
We  expect  a  significant  improvement  in  performance  if  true  Jacobians  are  used.  However,  in  that 
case  some  complications  may  arise  from  the  differentiation  of  the  Van  Albada  flux  limiter.  This 
matter  will  be  investigated  in  depth  in  future  research.  We  point  out  that  this  behavior  does  not 
contradict  the  ability  of  the  optimizer  to  handle  inexact  information.  The  optimizer  can  only  per¬ 
form  successfully  for  arbitrary  stopping  tolerances  if  the  degree  of  inexactness  can  be  adjusted  by 
the  optimizer  to  the  progress  it  makes  towards  computing  the  solution  and  thereby  to  the  required 
tolerance.  One  needs  to  adjust  the  stopping  tolerances  to  the  accuracy  in  function  values.  In  fact, 
we  could  have  relaxed  the  tolerance  on  the  coarse  grid  and  thereby  reduced  the  number  of  coarse 
grid  iterations,  while  maintaining  the  performance  on  the  finer  grids.  However,  since  an  exact  error 
bound  for  the  quality  of  the  Jacobians  used  is  not  available,  we  believe  that  one  tends  to  try  to  over¬ 
solve  the  problem.  Hence,  the  performance  displayed  in  Table  7.2  is  what  one  should  expect  in  the 
experimentation  phase  of  the  algorithm.  These  experiments  can  be  used  to  define  grid-dependent 
tolerances  which  will  lead  to  a  better  performance  than  that  shown  in  Table  7.2. 

We  used  the  optimizer  TRICE  because  of  its  capability  to  accept  solutions  of  linearized  state 
equations  and  adjoint  equations  computed  using  application  specific  solvers.  A  reformulation  of 
the  problem  presented  in  Section  6  was  necessary,  since  the  current  version  of  TRICE  only  solves 
problems  of  the  form  (5.1)-(5.3).  Some  inefficiencies  and  difficulties  might  be  attributed  to  this. 
It  is  expected  that  these  will  be  resolved  with  future  version  of  the  optimizer  for  solving  the  more 
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general  problems  (4.1)-(4.5).  A  high  percentage  of  the  number  of  block  matrix  solves  required 
to  solve  the  design  problem  (refer  Table  7.2)  can  be  attributed  to  the  large  number  of  iterations 
required  to  obtain  a  converged  solution  for  the  adjoint  equation.  This  is  especially  true  for  the  finer 
grids,  which  may  be  caused  by  ill-conditioning  in  the  grid.  We  return  to  this  issue  below.  The  high 
number  of  solves  can  also  be  partly  attributed  to  the  penalty  function  approach  we  use  to  address 
the  drag  constraint  (6.1).  When  the  drag  constraint  is  violated  the  gradient  of  the  objective  function 
with  respect  to  the  state  variables  becomes  very  large,  which  in  turn  means  that  the  residual  of  the 
adjoint  equation  is  very  large  and  requires  a  large  number  of  iterations  to  converge.  The  penalty 
term  (6.1)  also  causes  objective  function  to  change  abruptly  when  the  drag  constraint  is  violated. 
This  effect  is  currently  inadequately  reflected  in  the  model  (5.12)  used  to  compute  the  step  and  led 
to  a  large  number  of  unsuccessful  iterations. 

The  third  group  of  reasons  for  difficulties  in  the  solution  process  is  somewhat  related  to  the  first 
one  and  concerns  the  airfoil  design  problem  and  its  solution.  Jameson  [21]  has  shown  the  existence 
of  cases  where  nonunique  solutions  of  the  discretized  Euler  equations  can  be  obtained  for  certain 
airfoils.  While  the  all-at-once  approach  never  requires  the  solution  of  the  Euler  equation,  our  im¬ 
plementation  uses  the  solution  to  the  linearized  equations  and  adjoint  equations.  Existence  of  these 
solutions  and  their  dependence  upon  right  hand  side  data  need  to  be  investigated.  Another  possible 
reason  for  the  difficulties  in  convergence  behavior  is  the  fact  that  the  simple  algebraic  grid  that  we 
use  may  be  ill-conditioned.  Since  the  primary  purpose  of  the  present  study  was  to  demonstrate  the 
concept  of  using  the  all-at-once  approach  for  solving  the  design  problem,  we  have  not  investigated 
the  effect  of  the  grid  on  the  solution  process.  As  a  result,  no  attempt  has  been  made  to  ascertain  the 
quality  of  the  grid.  It  was,  in  fact,  observed  that  the  convergence  behavior  exhibited  by  the  flow 
solver  deteriorates  as  the  grid  is  refined,  which  could  be  an  indication  of  ill-conditioning.  Other 
grid  generation  techniques  and  airfoil  surface  discretizations  should  be  investigated  in  this  context. 
More  sophisticated  techniques  are  applied  in  many  of  the  papers  on  airfoil  design  cited  earlier  in 
this  paper.  However,  an  inclusion  of  such  techniques  in  the  all-at-once  approach  requires  a  care¬ 
ful  analysis  of  grid  sensitivities  which  are  needed  in  the  computation  of  Rw(q,w)  and  Jw{q.  w). 
Finally,  as  we  have  pointed  out  in  Section  2,  the  airfoil  design  problem  is  an  infinite  dimensional 
problem.  The  infinite  dimensional  problem,  its  discretization,  and  the  optimization  approach  have 
to  be  analyzed  jointly  to  derive  robust  and  efficient  solution  methods.  In  a  simplified  model  prob¬ 
lem  the  benefits  of  such  an  analysis  were  demonstrated  in  [6]  and  were  shown  to  lead  to  10-15% 
reductions  in  optimization  iterations.  Provisions  for  the  inclusion  of  infinite  dimensional  prob¬ 
lem  structure  into  the  optimizer  have  been  made,  see  [15].  Their  application  in  the  airfoil  design 
problem  are  part  of  future  research. 


8  Conclusion 

We  have  implemented  the  all-at-once  approach  to  solve  an  optimum  airfoil  design  problem.  The 
airfoil  design  problem  was  formulated  as  a  constrained  optimization  problem  in  which  flow  vari¬ 
ables  and  design  variables  are  viewed  as  independent  variables  and  in  which  the  coupling  steady 
state  2-D  Euler  equation  is  included  as  a  constraint.  To  implement  this  approach,  we  have  com- 
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bined  an  existing  optimization  algorithm,  TRICE,  with  an  existing  flow  code,  ErICA.  Details  of 
the  implementation  were  given  and  difficulties  arising  in  the  implementation  were  discussed.  Our 
numerical  results  indicate  that  the  cost  of  solving  the  design  problem  is  approximately  six  times 
the  cost  of  solving  a  single  analysis  problem.  This  is  consistent  with  the  expectation  that  the  de¬ 
coupling  of  flow  variables  and  design  variables  in  the  all-at-once  approach  makes  the  problem 
less  nonlinear  and  can  increase  the  efficiency  with  which  the  design  problem  is  solved.  Difficulties 
observed  in  the  solution  process  were  discussed  and  some  future  research  issues  were  addressed. 
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